Preparation

Packages

library(tidyverse) 
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(RColorBrewer)
library(treemapify)

Call custom function

source("Custom_function.R")

Read data

Read Nui Chua data with read.xlsx() function

#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# NuiChua plot
NCplot <- read.xlsx(filename, sheet = "NC_census2")

#Bidoup - NuiBa plot
#BDplot <- read.xlsx(filename, sheet = "BDNB_census2")

Display the Structure of data and view some first lines.

dbh: Diametre in mm;
lx,ly: tree coordinates to zero position of quadrat (20x20m);
gx,gy: tree coordinates to zero position of plot (100x100m).

str(NCplot)

## 'data.frame':    17383 obs. of  17 variables:
##  $ PlotID        : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ census        : chr  "C2" "C2" "C2" "C2" ...
##  $ tag           : chr  "01-0048" "01-0049" "01-0049" "01-0050" ...
##  $ stemtag       : chr  "01-0048_0" "01-0049_0" "01-0049_1" "01-0050_0" ...
##  $ spcode        : chr  "dimosp2" "crotth" "crotth" "dryppo" ...
##  $ quadrat       : chr  "0000" "0000" "0000" "0000" ...
##  $ lx            : num  1.119 1.545 1.545 1.188 0.327 ...
##  $ ly            : num  1.54 3.21 3.21 3.44 3.73 ...
##  $ gx            : num  1.119 1.545 1.545 1.188 0.327 ...
##  $ gy            : num  1.54 3.21 3.21 3.44 3.73 ...
##  $ dbh           : num  61 41 28 69 160 ...
##  $ codes         : chr  NA "TNBT" "TNBT" "TBV" ...
##  $ hom           : chr  NA NA NA NA ...
##  $ date          : chr  "2021-03-18" "2021-03-18" "2021-03-18" "2021-03-18" ...
##  $ family        : chr  "Euphorbiaceae" "Euphorbiaceae" "Euphorbiaceae" "Putranjivaceae" ...
##  $ genus         : chr  "Dimorphocalyx" "Croton" "Croton" "Drypetes" ...
##  $ scientificName: chr  "Dimorphocalyx sp.2" "Croton thorelii" "Croton thorelii" "Drypetes poilanei" ...

head(NCplot)

##   PlotID census     tag   stemtag  spcode quadrat       lx       ly       gx
## 1      1     C2 01-0048 01-0048_0 dimosp2    0000 1.118787 1.544663 1.118787
## 2      1     C2 01-0049 01-0049_0  crotth    0000 1.544927 3.205741 1.544927
## 3      1     C2 01-0049 01-0049_1  crotth    0000 1.544927 3.205741 1.544927
## 4      1     C2 01-0050 01-0050_0  dryppo    0000 1.188361 3.440553 1.188361
## 5      1     C2 01-0051 01-0051_0 syzysp2    0000 0.327383 3.727546 0.327383
## 6      1     C2 01-0052 01-0052_0 dimosp2    0000 1.127484 4.179776 1.127484
##         gy   dbh codes  hom       date         family         genus
## 1 1.544663  61.0  <NA> <NA> 2021-03-18  Euphorbiaceae Dimorphocalyx
## 2 3.205741  41.0  TNBT <NA> 2021-03-18  Euphorbiaceae        Croton
## 3 3.205741  28.0  TNBT <NA> 2021-03-18  Euphorbiaceae        Croton
## 4 3.440553  69.0   TBV <NA> 2021-03-18 Putranjivaceae      Drypetes
## 5 3.727546 160.5  <NA> <NA> 2021-03-18      Myrtaceae      Syzygium
## 6 4.179776 148.5  <NA> <NA> 2021-03-18  Euphorbiaceae Dimorphocalyx
##       scientificName
## 1 Dimorphocalyx sp.2
## 2    Croton thorelii
## 3    Croton thorelii
## 4  Drypetes poilanei
## 5      Syzygium sp.2
## 6 Dimorphocalyx sp.2

Community (Species composition, size-class structure and diversity)

Number of family, genus and species

Use n_distinct() counts the number of distinct value in column: family, genus, scientificName

NCplot |>  
  summarise("family" = n_distinct(family),
            "genus" = n_distinct(genus),
            "species" = n_distinct(scientificName))

##   family genus species
## 1     45    76      93

Number species by family

 (spByFam <- NCplot |> 
     group_by(family) |> 
     summarise(n = n_distinct(spcode)) |> 
     arrange(desc(n)))

## # A tibble: 45 × 2
##    family              n
##    <chr>           <int>
##  1 Euphorbiaceae       7
##  2 Rubiaceae           7
##  3 Phyllanthaceae      6
##  4 Annonaceae          5
##  5 Ebenaceae           5
##  6 Fabaceae            4
##  7 Lamiaceae           4
##  8 Celastraceae        3
##  9 Melastomataceae     3
## 10 Myrtaceae           3
## # ℹ 35 more rows

Figure in treemap layout

 library(treemapify)
ggplot(spByFam, aes(
  area = n,
  fill = n,
  label = paste(family, n, sep = "\n")
)) +
  geom_treemap(start = "bottomleft") +
  geom_treemap_text(place = "centre") +
  scale_fill_distiller(palette = "Greens", direction = 1) +
  labs(fill = "Species")

Species accumulation curves

Species accumulation curves is a graph recording the cumulative number of species recorded in a particular environment. It used to estimate the number of species in a particular area.

Here, we use b_specaccum() function.

curveNC <-  b_specaccum(NCplot)
head(curveNC)

##   subplot richness    S PlotName
## 1       1 35.24000  400         
## 2       2 48.11000  800         
## 3       3 56.08174 1200         
## 4       4 61.81360 1600         
## 5       5 66.23241 2000         
## 6       6 69.78524 2400

#FIGURE
ggplot(curveNC, aes(x = S, y = richness)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  geom_text(
    aes(label = round(richness, 1)),
    size = 4,
    vjust = -0.5,
    check_overlap = T
  ) +
  theme_bw(16) +
  labs(x = "Survey area (m2)", y = "Number of species")

Forest structure

We usually group trees into 3 size-class of DBH:

small tree [1-5)cm;
medium tree [5-20)cm;
big tree >=20 cm

NCplot <- NCplot |>
  mutate(
    dbh_cm = dbh / 10,
    # dbh value in mm
    gDBH = case_when(dbh_cm >= 20 ~ ">=20", 
                     dbh_cm >= 5 ~ "[5-20)", 
                     TRUE ~ "[1-5)")
  )

head(NCplot)

We count the number of stems at each diameter class.

(stemByDBH <- NCplot |>
    count(gDBH) |>
    mutate(freq = round(n / sum(n), 3)) |>
    arrange(desc(n)))

##     gDBH     n  freq
## 1  [1-5) 14221 0.818
## 2 [5-20)  2951 0.170
## 3   >=20   211 0.012

Figure

stemByDBH |>
  ggplot(aes(x = gDBH, y = n)) +
  geom_col(aes(fill = n)) +
  geom_text(aes(label = paste0(n, " (", percent(freq), ")")), vjust = -0.5) +
  scale_fill_viridis_c() +
  theme_bw(14) +
  labs(x = "DBH (cm)", y = "Nb of stem", fill = "Nb of stem")

We can also count the number of species at each diameter class.

(
  SpByDBH <- NCplot  |>
    group_by(gDBH) |>
    summarise(
      "family" = n_distinct(family),
      "genus" = n_distinct(genus),
      "species" = n_distinct(scientificName)
    )
)

## # A tibble: 3 × 4
##   gDBH   family genus species
##   <chr>   <int> <int>   <int>
## 1 >=20       26    33      39
## 2 [1-5)      41    71      87
## 3 [5-20)     38    59      72

Figure

SpByDBH |>
  tidyr::pivot_longer(-gDBH, names_to = "rank", values_to = "n") |>
  ggplot(aes(x = gDBH, y = n, fill = rank)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = n, label = n),
            position = position_dodge(width = 0.9),
            vjust = -0.5) +
  scale_fill_viridis_d() +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Rank")

Tree distribution map

Draw a blank plot with draw_blank_plot_1ha() function.

p <- draw_blank_plot_1ha()
p

Add add tree to map

p + geom_point(data = NCplot,
               aes(
                 x = gx,
                 y = gy,
                 size = dbh_cm,
                 color =  gDBH
               ),
               alpha = 0.7) +
  labs(color = "DBH class", size = "DBH(cm)")

For big tree

NCplotA <- NCplot |> filter(gDBH == ">=20")
p + geom_point(
  data = NCplotA,
  aes(x = gx, y = gy, size = dbh_cm),
  alpha = 0.7,
  shape = 21,
  fill = "blue",
  color = "white"
) +
  labs(size = "DBH(cm)")

For medium tree

NCplotB <- NCplot |> filter(gDBH == "[5-20)")
p + geom_point(
  data = NCplotB,
  aes(x = gx, y = gy, size = dbh_cm),
  alpha = 0.7,
  shape = 21,
  fill = "blue",
  color = "white"
) +
  labs(size = "DBH(cm)")

For small tree

NCplotC <- NCplot |> filter(gDBH == "[1-5)")
p + geom_point(
  data = NCplotC,
  aes(x = gx, y = gy, size = dbh_cm),
  alpha = 0.7,
  shape = 21,
  fill = "blue",
  color = "white"
) +
  labs(size = "DBH(cm)")

For a particular species, px: Dialium cochinchinense

## So do phan bo toan bo cay trong o mau
NCplotSP <- NCplot |> filter(scientificName == "Dialium cochinchinense")

p + geom_point(data = NCplotSP,
               aes(
                 x = gx,
                 y = gy,
                 size = dbh_cm,
                 color =  gDBH
               ),
               alpha = 0.7) +
  labs(color = "DBH class", size = "DBH(cm)")

We can also try with a detail DBH class

NCplot <- NCplot |> 
mutate(
  dbh_cm = dbh / 10 , # dbh in mm
  gDBH2 = case_when(
                dbh_cm >= 100 ~ ">=100",
                dbh_cm >= 90 ~ "[90-100)",
                dbh_cm >= 80 ~ "[80-90)",
                dbh_cm >= 70 ~ "[70-80)",
                dbh_cm >= 60 ~ "[60-70)",
                dbh_cm >= 50 ~ "[50-60)",
                dbh_cm >= 40 ~ "[40-50)",
                dbh_cm >= 30 ~ "[30-40)",
                dbh_cm >= 20 ~ "[20-30)",
                dbh_cm >= 10 ~ "[10-20)",
                TRUE ~ "[1-10)")
    )

stemByDBH2 <- NCplot |>
  count(gDBH2) |>
  mutate(freq = round(n / sum(n), 3)) |>
  arrange(desc(n))

stemByDBH2

##     gDBH2     n  freq
## 1  [1-10) 16210 0.933
## 2 [10-20)   962 0.055
## 3 [20-30)   162 0.009
## 4 [30-40)    33 0.002
## 5 [40-50)     6 0.000
## 6 [50-60)     6 0.000
## 7 [60-70)     2 0.000
## 8   >=100     1 0.000
## 9 [70-80)     1 0.000

stemByDBH2 |>
  ggplot(aes(x = gDBH2, y = n)) +
  geom_col() +
  geom_text(aes(label = n), vjust = -0.5) +
  #scale_fill_manual(values = c("grey","white"))+
  theme_bw(14) +
  labs(x = "DBHClass (cm)", y = "Number of stem")

Dominant species and Biodiversity Index

Dominant species

The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:

relative density of each species (RD)
relative basal area of each species (RBA)
relative frequency of each species (this can be omitted in 1 plot)

\[IVI={\left(relative.density + relative.basal.area \right)}\]

Basal area (BA) of each species can be calculated from diameter measured: $BA = pi * R^2$

Relative basal area (RBA): $RBA = BA / sum(BA)$

Relative density (RD) = number of individuals of each species / total

  NCplot <- NCplot |> 
  mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)

IVI <- NCplot |>
  group_by(scientificName) |>
  summarise(n = n(), sBA = sum(BA)) |>
  mutate(
    RBA = sBA / sum(sBA) * 100,
    RD = n / sum(n) * 100,
    IVI = RBA + RD,
    rank = rank(-IVI)
  ) |>
  arrange(desc(IVI))

Figure

ggplot(IVI, aes(x = rank, y = IVI)) +
  geom_line() +
  geom_point() +
  geom_label_repel(
    data = IVI[1:5, ],
    aes(label = paste0("italic('", scientificName, "')")),
    fill = NA,
    color = "black",
    size = 5,
    label.size = NA,
    min.segment.length = 0,
    segment.linetype = 2,
    segment.size = 0.4,
    parse = TRUE
  ) +
  theme_bw(14) +
  labs(x = "Species rank", y = "IVI (%)")

Diversity Index

A diversity index is a quantitative measure that reflects how many different types (such as species) there are in a dataset (a community). These indices are statistical representations of biodiversity in different aspects (richness, evenness, and dominance).

Richness

Richness is the total number of unique species

Shannon Diversity Index

As species richness and evenness increase, so diversity increases. The Shannon Diversity Index (H) is a way to measure the diversity of species in a community. The higher the value of H, the higher the diversity of species in a particular community. The lower the value of H, the lower the diversity. A value of H = 0 indicates a community that only has one species.

Evenness

The Shannon Equitability Index is a way to measure the evenness of species in a community. The term “evenness” simply refers to how similar the abundances of different species are in the community. This value ranges from 0 to 1 where 1 indicates complete evenness.

\[ Evenness = Shannon Index (H) / ln(richness) \]

Simpson’s Diversity Indices

Simpson’s Index (D) measures the probability that two individuals randomly selected from a sample will belong to the same species

With this index, 0 represents infinite diversity and 1, no diversity. That is, the bigger the value of D, the lower the diversity. This is neither intuitive nor logical, so to get over this problem, D is often subtracted from 1 to give:

Simpson’s Index of Diversity 1 - D

The value of this index also ranges between 0 and 1, but now, the greater the value, the greater the sample diversity. This makes more sense. In this case, the index represents the probability that two individuals randomly selected from a sample will belong to different species.

library(vegan)
# Simpson's Index of Diversity
  D <- diversity(IVI$n, index = "simpson")
# Shannon
  H <- diversity(IVI$n, index = "shannon")
# Richness
  S <- specnumber(IVI$n)
# Evenness
  E <- H/log(S)
# Create data frames 
  res <- data.frame("Richness"=S,"Simpson"=D,"Shannon"=H,"Evenness"=E)
  res

##   Richness   Simpson  Shannon Evenness
## 1       93 0.9014608 2.927733 0.645928

For each DBH class

Create function

Create calcIVI() to calculator Importance Value Index and figureIVI() to for plotting figure

calcIVI <- function(dataplot) {
  IVI <- dataplot |>
    mutate(G = ((dbh / 2) / 1000)^2 * pi) |>
    group_by(scientificName) |>
    summarise(n = n(), sG = sum(G)) |>
    mutate(
      UT = sG / sum(sG) * 100,
      MD = n / sum(n) * 100,
      IVI = UT + MD,
      rank = rank(-IVI)
    ) |>
    arrange(desc(IVI))
  return(IVI)
}

figureIVI <-  function(IVI) {
  ggplot(IVI, aes(x = rank, y = IVI)) +
    geom_line() +
    geom_point() +
    geom_label_repel(
      data = IVI[1:5, ],
      aes(label = paste0("italic('", scientificName, "')")),
      fill = NA,
      color = "black",
      size = 5,
      label.size = NA,
      min.segment.length = 0,
      segment.linetype = 2,
      segment.size = 0.4,
      parse = TRUE
    ) +
    theme_bw(14) +
    labs(x = "Species rank", y = "IVI (%)")
}

Create biodivIndex() to calculator Diversity Index

biodivIndex <- function(IVI) {
  H <- diversity(IVI$n, index = "shannon")
  D <- diversity(IVI$n, index = "simpson")
  S <- specnumber(IVI$n)
  E <- H / log(S)
  res <- data.frame(
    "Richness" = S,
    "Simpson" = D,
    "Shannon" = H,
    "Evenness" = E
  )
  return(res)
}

Run

Big tree (DBH >= 20 cm)

IVIA <- calcIVI(NCplotA)
head(IVIA)

## # A tibble: 6 × 7
##   scientificName             n    sG    UT    MD   IVI  rank
##   <chr>                  <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Syzygium sp.1             24 3.20  22.8  11.4   34.1     1
## 2 Syzygium sp.2             37 2.05  14.6  17.5   32.1     2
## 3 Drypetes poilanei         24 1.05   7.48 11.4   18.9     3
## 4 Diospyros sp.2            21 1.19   8.44  9.95  18.4     4
## 5 Dialium cochinchinense    18 1.11   7.87  8.53  16.4     5
## 6 Dehaasia sp.              10 0.753  5.36  4.74  10.1     6

figureIVI(IVIA)

biodivIndex(IVIA)

##   Richness   Simpson  Shannon  Evenness
## 1       39 0.9170055 2.919805 0.7969853

Medium tree (DBH in 5-20 cm)

IVIB <- calcIVI(NCplotB)
head(IVIB)

## # A tibble: 6 × 7
##   scientificName        n    sG    UT    MD   IVI  rank
##   <chr>             <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Syzygium sp.2       549 5.64  26.0  18.6  44.6      1
## 2 Drypetes poilanei   521 3.76  17.3  17.7  35.0      2
## 3 Diospyros sp.2      227 1.83   8.43  7.69 16.1      3
## 4 Cleistanthus sp.    327 1.08   4.97 11.1  16.1      4
## 5 Walsura bonii       136 0.868  4.00  4.61  8.61     5
## 6 Diospyros sp.1      116 0.853  3.93  3.93  7.86     6

figureIVI(IVIB)

biodivIndex(IVIB)

##   Richness   Simpson  Shannon  Evenness
## 1       72 0.9062302 2.945215 0.6886708

Small tree (DBH in 1-5 cm)

IVIC <- calcIVI(NCplotC)
head(IVIC)

## # A tibble: 6 × 7
##   scientificName            n    sG    UT    MD   IVI  rank
##   <chr>                 <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Cleistanthus sp.       2767 1.52  22.8  19.5  42.2      1
## 2 Dimorphocalyx sp.2     3215 1.10  16.5  22.6  39.1      2
## 3 Goniothalamus sp.      2230 1.13  17.0  15.7  32.6      3
## 4 Drypetes poilanei       599 0.465  6.97  4.21 11.2      4
## 5 Cleistanthus monoicus   555 0.340  5.09  3.90  8.99     5
## 6 Mallotus sp.            726 0.176  2.63  5.11  7.74     6

figureIVI(IVIC)

biodivIndex(IVIC)

##   Richness   Simpson  Shannon  Evenness
## 1       87 0.8767146 2.706991 0.6061457

Biomass and carbon stock

Biomass

Above-ground biomass (AGB)
Below-ground biomass (BGB)

AGB calculation

Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87–99.

AGB=Total aboveground biomass (kg)
D=DBH (cm)
H=Height (m)
ρ=Wood density (g/cm3)
BA=Basal area (cm2)

	Equation
Tropical forest DRY	$AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)$
	$AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$
Tropical forest MOIST	$AGB = exp(-2.977 + ln(ρ * D^2 * H))$
	$AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$

library(BIOMASS)

## For more information on using BIOMASS, visit https://umr-amap.github.io/BIOMASS

## Using temporary cache
##   It is recommended to use a permanent cache to avoid to re-download files on each session.
##   See function createCache() or BIOMASS.cache option.

woodensity <- getWoodDensity(
  genus = NCplot$genus,
  species = NCplot$scientificName,
  family = NCplot$family
)

## The reference dataset contains 16467 wood density values

## Your taxonomic table contains 93 taxa

head(woodensity)

NCplot$meanWD <- woodensity$meanWD

$AGB=ρ∗exp(−0.667+1.784∗ln(D)+0.207∗(ln(D))2–0.0281∗(ln(D))3)$

NCplot <- NCplot %>% 
  mutate(dbh_cm = dbh/10,
         AGB = meanWD * exp(-0.667 + 1.784*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3))

(AGB <- sum(NCplot$AGB)/1000) #Mg/ha

## [1] 264.7534

BGB calculation

The ratio developed by Mokany et al. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level

(BGB = 0.275 * AGB)

## [1] 72.80718

Total biomass of 1 hecta

(TB <- AGB + BGB)

## [1] 337.5605

Carbon stock

The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.

#Carbon stock (ton)
(TC <- TB * 0.47)

## [1] 158.6535

# CO2 (ton)
(CO2 <- TC * 44/22)

## [1] 317.3069

Carbon credit

1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.

In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)

1 hecta forest at Núi Chúa national park ~ 1586 $

(CO2 * 5)

## [1] 1586.535

	Equation
Tropical forest DRY	\(AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)\)
	\(AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\)
Tropical forest MOIST	\(AGB = exp(-2.977 + ln(ρ * D^2 * H))\)
	\(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\)

Day 01 - Nui Chua Forest Dynamics Plots

Nguyen-Le Xuan Bach

2025-11-30

Preparation

Packages

Call custom function

Read data

Community (Species composition, size-class structure and diversity)

Number of family, genus and species

Number species by family

Species accumulation curves

Forest structure

Tree distribution map

Dominant species and Biodiversity Index

Dominant species

Diversity Index

Richness

Shannon Diversity Index

Evenness

Simpson’s Diversity Indices

For each DBH class

Create function

Run

Biomass and carbon stock

Biomass

AGB calculation

BGB calculation

Carbon stock

Carbon credit